
#include <complex.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "f2c.h" 
#include <sys/time.h>
#include <malloc_align.h>

#include "cmatrixprod.h"
#include "cmatrixprod_ref.h"

#define PERROR(s) {fprintf(stderr,s);exit(1);}

//#define AFFICHE 1


void affiche(t_complex a[], int nlig, int ncol, char *s)
{
  int i,j;
  printf("\n\n%s=[\n",s);
  
  for (i=0; i<nlig;i++) {
    for (j=0;j<ncol;j++) 
      printf("%f + %f*i  ",a[i*ncol+j].r,a[i*ncol+j].i);
    printf(";\n"); 
  }   
 printf("];\n");
}


float myrand()
{
  return(floor(rand()/(float)RAND_MAX*10)+1);
}


double my_gettimeofday()
{
  struct timeval tmp_time;
  gettimeofday(&tmp_time, NULL);
  return tmp_time.tv_sec + (tmp_time.tv_usec * 1.0e-6L);
}



int main(int argc, char*argv[])
{
  t_complex *a, *b, *c, *cc;
  int i,j;
  int m,n,p;
  int erreur=0;
  double t1, t2;

  if (argc!=4)  PERROR("testzgemm : m, n, p (matrice a[m][p], b[p][n]");

  m=atoi(argv[1]);
  n=atoi(argv[2]);
  p=atoi(argv[3]);
  
  printf("m %d    n %d      p %d\n",m,n,p);

  srand(16);
  
  if ((a=(t_complex *)_malloc_align(sizeof(t_complex)*m*p,4))==NULL)
    perror("allocation a\n");

  if ((b=(t_complex *)_malloc_align(sizeof(t_complex)*p*n,4))==NULL)
    perror("allocation b\n");

  if ((c=(t_complex *)_malloc_align(sizeof(t_complex)*m*n,4))==NULL)
    perror("allocation c\n");

  if ((cc=(t_complex *)_malloc_align(sizeof(t_complex)*m*n,4))==NULL)
    perror("allocation cc\n");


  for (i=0; i<m*p;i++) {
    a[i].r=myrand();
    a[i].i=myrand();
  }
  //#ifdef AFFICHE 
  affiche(a,m,p,"a");
  //#endif


  for (i=0; i<p*n;i++) {
    b[i].r=myrand();
    b[i].i=myrand();
  }
#ifdef AFFICHE 
  affiche(b,p,n,"b");
#endif



  for (i=0; i<m*n;i++) {
    c[i].r=1;
    c[i].i=i%2;
    cc[i].r=1;
    cc[i].i=i%2;
  }


  
  t1=my_gettimeofday();
  cmatrixprod(m, n, p, 	a, b, c);
  t2=my_gettimeofday();

  
  cmatrixprod_ref(m, n, p, 	a, b, cc);


#ifdef AFFICHE 
  affiche(c,m,n,"c");
#endif

#ifdef AFFICHE 
  affiche(cc,m,n,"cc");
#endif


  printf("Resultat:\n");

  for (i=0; i<m;i++){
    for (j=0;j<n;j++){
      if (fabsf((c[i*n+j].r-cc[i*n+j].r)/cc[i*n+j].r) >1.e-5  ||
	  fabsf((c[i*n+j].i-cc[i*n+j].i)/cc[i*n+j].i) >1.e-5 ){
	erreur=1;
	printf("Probleme ligne %4d colonne %4d  calculé : %.7e + %.7e i  attendu calculé : %.7e + %.7e i  \n",
	       i,j,c[i*n+j].r,c[i*n+j].i, cc[i*n+j].r,cc[i*n+j].i);      
      }
    }
  } 
  if (!erreur){
    printf("PROGRAMME CORRECTE\n");
  }else{
    printf("Programme incorrecte\n");
  }

  printf("Temps en seconde de votre fonction  %.3lf\n",t2-t1);
  printf("Puissance de calcul en GFLOPS = %.6lf\n",8.*m*n*p/(t2-t1)/1.e9);

  return 0;

}
  

    


